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Abstract. We study photon condensation phenomena in a driven and dissipative 
array of superconducting microwave resonators. Specifically, we show that by using 
an appropriately designed coupling of microwave photons to superconducting qubits, 
an effective dissipative mechanism can be engineered, which scatters photons towards 
low- momentum states while conserving their number. This mimics a tunable coupling 
of bosons to a low temperature bath, and leads to the formation of a stationary photon 
condensate in the presence of losses and under continuous-driving conditions. Here we 
propose a realistic experimental setup to observe this effect in two or multiple coupled 
cavities, and study the characteristics of such an out-of-equilibrium condensate, which 
arise from the competition between pumping and dissipation processes. 
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1. Introduction 

Massive particles and spin systems have been at the center of investigations in quantum 
many-body physics since the inception of the field. Photonic systems, however, have 
been the object of interest in this field only much more recently. This is, on the one 
hand, due to the fact that photons are intrinsically non-interacting, and on the other 
hand, due to the difficulty to confine and experiment with photons before they decay. 
One important exception in this context are ideas to produce a Bose-Einstein condensate 
(BEC) of weakly interacting polaritons [H EJ [3J [H [5], whose large photonic component 
is essential to achieve a low effective mass and therefore transition temperatures which 
are accessible in solid-state experiments. While current cold-atom experiments offer a 
well controlled setting to study bosonic quantum gases [6] , the investigation of photonic 
condensates is still an active area of research, stimulated by new experimental directions 
in the field [H El [9] and controversial discussions on the exact nature of such condensates 

m HD H2J- 

More recently, the connection between many-body physics and photons has been 
addressed again from the perspective of photonic quantum simulators [HI [15l [HI 
[T71 [T8l [T9l 1201 I2T1 122] . Here, the general goal is to implement strongly interacting 
many-body systems in a controllable way, to simulate the properties of non-trivial 
condensed-matter models. This can be achieved, for example, using ideas from cavity 
quantum electrodynamics (QED) [23], where effective photon-photon interactions can be 
obtained through the coupling to an intermediary system, such as atoms [211 [25] or solid- 
state systems [26], [271 EH]. Based on this principle, various schemes for implementing 
bosonic Hubbard models for photons on a lattice pH EL photonic quantum Hall 
systems [29l EQl EH [32], and strongly interacting photons in a ID continuum j33j [34] . 
have been theoretically investigated. Although the experimental implementation of 
these ideas is still challenging, the development of scalable cavity-QED systems in on- 
chip devices [35] is rapidly progressing, and analogous ('circuit QED') systems in the 
microwave regime j36j EH EEJ EHJ HDJ [41] already approach the stage where photonic 
lattice models can be realized. 

A central and still open question in the field of photonic quantum simulators is 
how to prepare and probe quantum many-body states of light. Because of unavoidable 
losses and the absence of a chemical potential (arising from an equilibrium particle 
reservoir), photonic many-body systems must necessarily be studied under driven 
and non-equilibrium conditions. Therefore, familiar concepts from condensed-matter 
physics, such as ground-state or equilibrium phase diagrams, are a priori not accessible 
in photonic many-body systems. In previous works it has been suggested to study, 
for example, the transient dynamics of an initially pumped system [T51 CEH [42], or to 
use excitations spectroscopy and photon-correlation measurements of a weakly driven 
system to infer certain properties of the interacting photonic system [HU [43], EH 1451 146] . 
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Figure 1. Schematics of a coupled array of photonic cavities (represented in grey). 
The photons in each cavity are tunnel-coupled to neighboring sites with amplitude 
J and decay with rate T. The cavities are driven with an external coherent field of 
strength This external driving compensates losses and ensures a finite stationary 
photon population. 



In this work we consider a different scenario and study the dissipative dynamics 
of a strongly driven cavity array in the presence of engineered dissipation [471 SHI 
1491 l50l I5T1 1521 1531 154] . More precisely, we will show how in a circuit QED setting, 
the coupling of microwave photons to superconducting qubits can be used to scatter 
photons from high to low momentum states, eventually accumulating photons in a 
zero-momentum condensate. Previously, a similar scheme has been explored as a 
dissipative way to prepare a Bose Einstein condensate (BEC) of atoms [48]. In the 
case of photonics systems, it is essential that this mechanism implements a number- 
conserving coupling of photons to an effective low-temperature bath, and therefore 
provides a new approach for the preparation of quasi-equilibrium many-body states in 
open and driven photonic systems. Here we propose and analyze a proof-of-principle 
experiment to study condensation of microwave photons in a dissipative cavity array, and 
describe the properties of such a non-equilibrium BEC, which arise from the interplay 
between driving, decay, and thermalization processes. Such experiment would realize a 
controlled setting for the simulation of non-equilibrium condensation phenomena, and 
more generally, would offer a new route towards preparing stationary states of strongly 
interacting photons. Moreover, our work shows that the flexibility provided by circuit 
QED, which so far has mainly been employed to engineer strong coherent interactions, 
can equally well be used for the design of various non-trivial dissipative couplings, and 
be applied to simulate driven quantum systems in unconventional environments. 

2. Coupled cavity arrays 

Fig. [1] illustrates the basic setup for an array of L coupled cavities, where each cavity is 
represented by a single photonic mode of frequency u c and bosonic annihilation operator 
C£. Photons can tunnel between neighboring cavities with hopping amplitude J and, in 
addition, local interactions with two- or many-level systems can be used to induce Kerr- 
type nonlinearities with an effective photon-photon interaction strength U. A generic 
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model for this system is then given by the Bose-Hubbard Hamiltonian (see e.g. Ref. [13J) 

H c = Y]cj c c\ci - J J^(qcJ+i + Q+icJ) + ^ ^2 AAwi- (!) 

tt t 

Various generalizations of this model have been discussed in the literature and several 
implementations have been proposed using optical or microwave cavities. The photonic 
Hubbard model given in Eq. (OQ) describes a gas of interacting bosons on a lattice, where 
the hopping J competes with on-site interaction U. However, in contrast to the cold- 
atom physics scenario, the bosons here can decay, and under experimentally-relevant 
conditions uj c 3> J, [/, k^T (where T is the temperature and &b the Boltzmann constant), 
the equilibrium state of this model is simply the vacuum state. Therefore, in photonic 
many-body systems we are mainly interested in the out-of-equilibrium dynamics of % c 
in the presence of losses and external driving fields. In particular, in this work we 
model the resulting dissipative dynamics for the system density operator p by a master 
equation (ME) of the form, 

p = -i[H c + Hn(t),p] + T ^ V[c £ ]p + C K p, (2) 

£ 

where V[c]p = 2cp& — &cp — p&c. In Eq. ([2j) the Hamiltonian T~Ln(t) = J2 £ Qt(e~' luJdt c\ + 
e iUJdt Ci) describes an external driving field of frequency which is used to excite the 
system, and the second term accounts for photon losses in each cavity with a field decay 
rate Y. While a finite driving field is required to counteract the losses, it will in general 
also compete with % c and, for strong driving fields, even dominate the system dynamics. 
Therefore, in previous works it has been suggested to study either the transient dynamics 
of an initially prepared photonic state p~5j [EH [42] (where = for times t > 0) [42] or 
to use weak excitation spectroscopy [311 S31 HH SSJ [46] (f|g 0) to probe the many-body 
spectrum of Hamiltonian H c . 

In this work we are interested in the opposite regime of a strongly and continuously 
driven cavity array. We study the dynamics of this system in the presence of an 
additional artificial thermalization mechanism, denoted by C K in Eq. ([2]). More precisely, 
we will show below how a non-local coupling of photons to superconducting qubits can be 
engineered in an array of microwave cavities to implement a dissipative photon scattering 
process of the form 

= E + a 5+i)(^ - M] + - + q + i)]. (3) 

£ 

The interpretation of this term can be seen best in the case of just two cavities. Then, 
the first term in Eq. ([3]) describes the scattering of photons from the asymmetric 
(energetically higher) mode c a = (ci — c 2 )/V^ into the symmetric (energetically lower) 
mode c s = {c\ + C2)/v / 2, while preserving the total photon number. The second term in 
Eq. (|3]) describes the reverse process. In the absence of losses, the action of C K would, for 
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two cavities, drive arbitrary initial photon states into a thermal equilibrium distribution, 
characterized by a detailed balance between symmetric and antisymmetric modes, with 
the ratio k 1 /k = exp(— 2J//cBT e ff) defining an effective temperature T e ff. Similarly, for 
a whole cavity array, the processes ~ k scatter photons to lower momenta and - as 
previously shown for an analogous system of cold atoms [48] - drive the system towards 
a condensate of photons in the zero momentum mode, |BEC) = (Cq =o ) N \0) / v^VT, where q 
is the quasimomentum and TV the total number of particles. Also in this case the opposite 
processes ~ k! can be roughly interpreted as a finite temperature effect, although, as it 
will be shown below, for L > 2 and k! ^ the stationary state of C K does no longer 
represent a true thermal equilibrium state. 

The competition of C K with a Kerr-type nonlinearity in the Hamiltonian has already 
been discussed in the atomic, number-conserving setup [55l [56] . However, in the 
photonic case, under realistic conditions, C K competes with external driving fields and 
intrinsic photon losses as described by the full ME ([2j) . To study this competition more 
explicitly, we focus in the following on the limit U 0, and consider the experimental 
scenario where neighboring cavities are driven by coherent fields of alternating amplitude 
Qi = (—iyQ (see Fig. Q]). In this case, the accumulation of photons in the q = 
mode can be interpreted as a clear signature for photon condensation induced by the 
dissipative photon-photon interactions in C K . 

3. Physical implementation 

In this section we show how the photon scattering mechanism described by Eq. ([3]) 
can be implemented in an array of superconducting microwave resonators. The basic 
setting is illustrated in Fig.EK for two cavities, each coupled to a nonlinear element, for 
example a superconducting Cooper pair box ('charge qubit'). The qubits are placed next 
to each other to obtain strong electrostatic or magnetic interactions. As we discuss now, 
this configuration allows us to engineer both nonlinear as well as non-local dissipation 
processes for photons. 

For simplicity we restrict the following analysis to a single block of only two cavities, 
but the generalization to a whole array of linked cavities is straightforward. The 
Hamiltonian for this system can be derived from the corresponding equivalent circuit 
model schematically shown in Fig. [2^. By restricting each resonator to a single mode 
q, £ = 1,2, and by approximating each Cooper pair box by a two-level system with 
ground state \g) and excited state |e), we obtain 

w=£ "4* + J2 a4*- + ) + n q (t), (4) 

£=1,2 £=1,2 

where = The first part is the free cavity Hamiltonian, while the second term 

describes the dipole interaction between photons and qubits with strength g. Finally, 
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Figure 2. Dissipation engineering with circuit QED. a) Implementation of the 
Hamiltonian in Eq. (j4j). Two stripline resonators are coupled through a system of 
two mutually interacting Cooper pair boxes, b) Four- level diagram corresponding to 
the two coupled qubit s, whose dynamics is described by Eq. (|6j). 



7i q (i) is the Hamiltonian of the two coupled qubit s, which we assume of the generic 
formQ 

£=1,2 

+ XA 1 Mf + K(a^ + a^). (5) 

Here, Q<p are the Rabi frequencies of an applied microwave field, which is used to drive 
the qubits at frequency ujq. The coupling constants X z and X x are the strengths of the 
diagonal and off-diagonal qubit-qubit interactions, respectively. The Hamiltonian (|4j) 
can be rewritten in terms of the symmetric c s = [c\ + C2)/v / 2 and antisymmetric 
c a = (ci — 02) /y/2 cavity modes, and the qubit states \E) = |ee), \G) = \gg), 
\S) = (\eg) + \ge))/y/2, and \A) = (\eg) — \ge))/y/2. In this new basis, choosing 
fiq 1 ^ = — f^q 2 ^ = Q a /y/2 : and changing into a frame rotating with cj , we obtain 

H = - 5c% - 8c\c a - A e \E)(E\ - A S \S)(S\ - A a \A){A\ 
+ n a [\A){G\-\E)(A\ + R.c] 

+ g[ct(\G)(S\ + \S)(E\) + ci(\G)(A\-\A)(E\) + K.c.}, (6) 

where 5 = uj — u C) A s?a = uj — uj q =p X x , A e = 2(uj — uj q ) — X z . Apart from the coherent 
dynamics described by Eq. ([6j) , we include dissipation in the form of intrinsic cavity loss 

f The explicit derivation using circuit theory is omitted here. For a comprehensive review on the topic 

c.f. [SSlEflEBl. 
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with rate 2Y and qubit decay with rate 7. The full system dynamics is then described 
by the ME 

p = -i[u, P ]+ £ (rvic^p + ^viwi^p + ^vuiEttp). (7) 

7]= 3,0, 

and a summary of the energy levels and the most relevant transitions is shown in Fig.Eb. 

Our goal in the following is to eliminate the qubit dynamics and derive an effective 
ME for the cavity modes. Specifically, we are interested in the dissipative two-photon 
process, where the qubits change from \A) to \S) by absorbing a photon from the 
antisymmetric cavity mode (c a ) and emitting a photon into the symmetric cavity mode 
(cj) (see Fig. Eb). Since the overall process is 0(g 4 ) in the photon-qubit coupling 
strength, a general derivation is quite involved, and to make the following discussion 
more transparent we will now focus explicitly on the hierarchy of energy scales drawn 
in Fig. Eh- In particular, we assume that 5, A e , and A = 5 — A e are much larger than 
all the other frequency scales g, Q a , A s?a , T, 7. In this limit, none of the single-photon 
processes is resonant and we can use a Schrieffer- Wolff transformation to derive an 
effective two-photon Hamiltonian. We perform the unitary transformation % = VHV^ 
with V — e s : and make the ansatz 

S = ct [a ltS \G)(S\+a 2tS \S)(E\\ +c\ [a 1>a \G) {A\ - a 2 , a \A) (E\\ -H.c. (8) 

We define and % g as the first and last lines of the Hamiltonian ([6]), respectively, and 
choose the coefficients such that Hq] = —'Hg. This can be achieved by setting 

In view of the assumptions discussed above, we will use the approximate results 
^1,77 ~ g/d, a 2,r] ~ 9 1 'A and define the parameter e = A/ 5 < 1. After this transformation 
the Hamiltonian (|6j) reads 

H = H q + H c + H int + 0{g\ gQ a ). (10) 

Here "H q is the qubit Hamiltonian 'H q written above, with small frequency shifts absorbed 
into a redefinition of the detunings A e s a . The modified cavity Hamiltonian is 

H c = -8c\c s -8c\c a -J(c\c s -clc a )+g eS [{c\c a + ec\c s ){P sa ) + H.c] , (11) 
where g etf = g 2 /A, 5 = S - g etf ({P aa ) + (P M »/2, J = g eS ({P aa ) - (P M »/2, P sa = \S){A\, 
P ss = (\S)(S\ - \E)(E\) + e(\G)(G\ - \S)(S\), (12) 
Paa = {\A){A\ - \E)(E\) + e(\G)(G\ - \A){A\), (13) 

and the average (•) is taken with respect to the stationary qubit state in the limit 
g —> 0. For the parameter regime of interest, (P sa ) ~ 0, and % c corresponds to the free 
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cavity Hamiltonian with a qubit-mediated tunneling amplitude J. Finally, the effective 
coupling between photons and qubits can be written as 

Hint - 9e$ [clcsPss + c\c a P aa + {c% + €c\c s )P sa + {c% + ec\c s ) ] P^] , (14) 

where Pk m = Pkm — (Pkm)- Note that in 7i int we have only kept resonant two-photon 
processes and already omitted terms like ~ gtt a / A x c 5?a and g 2 / A x c 2 sa . Both of these 
terms oscillate with the detuning S and will only give small corrections to the results 
presented below. 

The transformation V induces a weak mixing between photon and qubit states, 
which apart from generating new effective interactions also modifies the dissipative 
couplings. Since we consider 7 > T, we find that the main result of this mixing is an 
enhancement of the cavity decay rates T —> T a s = T + f s?a . After a transformation of 
the jump operators in the ME (ED, and averaging over the qubit states, we find that for 
the parameter regime of interest, these rates are 

t~^ 2 e 2 [(\G)(G\) + {\A)(A\)}, (15) 

f « ~ ^ b 2 (\G)(G\) + (2 + e 2 )(|^)(^|)] . (16) 

In summary, we find that the system dynamics in the new dressed state 
representation is well described by the ME 

P=(£ q + £ c + C- mt )p : (17) 

where 

C qP = -t[H q , p] + J2 (lnG){v\]p+lnv)(E\]p), (18) 

r]=s 1 a 

C cP = -i[H c , p]+J2 r v v i d r,}p, (19) 

r]=s,a 

£intP = -i[H [nU p}. (20) 

In the limit g e R — >■ 0, qubits and cavities are decoupled, and the system relaxes on a 
timescale 7 -1 into the state p(t) ~ p c (t) ® where £ q p^ = 0. For finite g e $ < 7, we 
can use a perturbative projection operator technique to eliminate the qubit degrees of 
freedom and derive an effective ME for the reduced cavity density operator [59], 

Pc = £cP- J drTT q {[H int) e rC - ([U intl p c (t) <g> pj])]}. (21) 

Evaluating this expression we obtain 

p c « - i[U c + U e ft, p c ] + ^2 T v V[c v ]p c + nV[(c\c a + ec\c s )]p c 

r]=s 1 a 

+ ^2 r W i^pcfir,' + h v ,p c h v - h v h v ,p c - p c n v n v >) , (22) 
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where = c^c^ and is an effective photon-photon interaction 

Heff = £4ff(C s 4 + ^ a c\){c\c a + 6c\c s ) + ^ U^fl^ . (23) 

r),r}'=s,a 

In Eq. ([22D we have defined = ^ ff Re{5 aS5Sa (0)}, T^, = 0eff Re {£w7'*7'(°)}> and the 
interactions C/ e ff = g^Iml^^O)} and , = g^ImlS^^/^O)}, given in terms of the 
qubit correlation spectra 

roc 

Srf rfre-(^(r)P«,(0)), (24) 

Jo 

which can be calculated using the quantum regression theorem [50] . 
3.1. Discussion 

By ignoring coherent photon-photon interactions for the moment, we see that the first 
line of the effective cavity dynamics given in Eq. f |22l) represents the dissipatively coupled 
cavities as introduced in Eq. ([2j). In particular, we can write 

kV[{c% + ec\c s )\ = KV[c\c a ] + K(?V[c\c 8 ] + C e . (25) 

By comparison of Eqs. f l22l ) and ([3]) we identify k! = e 2 ft, with the parameter e 2 = (A/5) 2 
determining the minimal effective temperature of the engineered two-photon process. 
For typical parameters, u q = 10 GHz, S = 1 GHz and A e = 0.8 GHz, we obtain e = 0.2, 
and hence the limit k! considered in most parts of the paper can be implemented 
to a very good approximation. However, we point out that the ratio k' / k can always be 
increased artificially, for example, by adding an additional coherent or incoherent driving 
field to populate the state \S). Our derivation also shows the existence of coherences 
between c\c a and c\c S) which lead to additional squeezing terms of order 0{ne) ) and 
that in Eq. f |25l) are summarized by C e . For e«l, and in the presence of cavity losses, 
we expect the influence of these coherences to be negligible, while the relevant effect on 
the populations is already captured by a finite n f . Indeed, numerical simulations similar 
to those discussed in Sec. [4] show no significant differences between the exact ME f |22l) 
and our model Liouvillian (|3]), even up to values of e « 0.5. However, by choosing e « 1 
and T <C ft, it would be possible to extend our model and study physical effects directly 
linked to the parameter e, such as number-conserving squeezing as previously analyzed 
in [53. 

Apart from the two-photon scattering processes of interest, the coupling to the 
qubits introduces various imperfections. On the one hand, these are the enhanced 
cavity decay rates T s a introduced above, and on the other hand, we obtain additional 
dephasing terms Tfj due to fluctuating Stark shifts when the qubit jumps incoherently 
between states \G) and \A). Since the ratio n/Y s ^ a scales as ^ 5 ,2 /7 2 , single photon 
losses can be suppressed in the strong coupling regime, as long as ft also exceeds the 
bare decay rate Y. The ratio i^/Tf- is independent of g, and depends on the details 
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Figure 3. a) Single and two-photon rates which appear in the effective model of 
Eq. ([22]) . The photon scattering rate of interest, k, becomes significantly larger than 
the others in a certain region of the parameter space. Here we have chosen g = 25 MHz, 
S = 1 GHz, A s = A a = MHz, A e = 800 MHz, 7 = 10 MHz, r = 10 kHz. b) Figure 
of merit for the validity of our model for high-momentum modes, F a = K,/(T a + rj a ), 
as a function of qubit driving Q a and detuning A a . Rest of parameters as in (a), c) 
The same plot as in (b) but for low-momentum modes, F s = k/(T 3 + rf a ). 



of the correlation functions f |24l ) and therefore on the parameters Q a , A s?a . In Fig. EH, 
ft is compared with the other decoherence rates, and we see that an experimentally 
feasible parameter range exists, where k is the dominant rate. The main corrections 
arise from the enhanced decay T a and the dephasing T^ a of the antisymmetric mode, 
while the corresponding rates for the symmetric mode are much smaller. In Fig. Eb and 
Eb we plot the two quantities F v = n/iT^ + I^L), r\ = s, a, as a function of the control 
parameters f2 a , A a . These quantities represent a simple figure of merit for characterizing 
the validity of our model for low (F s ) and high (F a ) momentum modes. We see that 
within a large parameter regime, the number conserving photon scattering rate k can 
exceed all decoherence processes discussed here, and therefore the proposed model in 
Eq. (Ej) provides a reliable description of the system dynamics. In particular, this is true 
for the low momentum (symmetric) modes, where condensation occurs. 

Finally, let us briefly comment on the effective photon-photon interactions described 
by Heff- For a resonant two-photon process, A s = 0, we find U e $ « 0. Moreover, 
U® v , w , < k for the typical parameters considered above, and in particular these 
interactions vanish for A a = 0. Therefore, in this work we focus exclusively on the 
purely dissipative cavity dynamics. However, by setting |A 5 | > 7 instead, the regime of 
strong coherent photon-photon coupling c7 eff > k > I^L, T s a can be engineered as well, 
having ~ U&c\c s c\c a . 
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In summary, we find that the ME ([2j) can be implemented with superconducting 
microwave cavities under realistic experimental conditions. Our analysis, which has 
been detailed here for generic coupled two level systems, can be thus adapted to various 
nonlinear Josephson devices, where the parameter X XjZ plays the role of the nonlinearity. 

4. Condensation of photons: two cavities 

4-1- Semiclassical approximation 

In the previous section we have shown how to implement the effective ME (Ej). Now, 
we discuss how the scattering between photons induced by the Liouvillian C K leads to 
condensation of photons in an open and driven cavity array. In this section we first 
illustrate the basic process for the minimal case with L = 2 cavities only. As already 
mentioned above, we assume [7 = 0, and consider the case where only the antisymmetric 
mode is excited, i.e. Qi = —Vt 2 = In the rotating frame with respect to the drive 

frequency uj^ the equations of motion (EOM) for the field expectation values (c s ) and 
(c a ) of the symmetric and antisymmetric cavity modes are 



where 5^ = uj c + J — uj& is the detuning of the antisymmetric mode with respect to the 
driving field. For the moment, we assume k! = in Eq. ([3]) and postpone the discussion 
of this term to the end of the present section. The EOM for the field expectation 
values couple to third-order correlation functions, starting a hierarchy of equations that 
cannot be solved analytically. To truncate the hierarchy, we resort here to a semiclassical 
approximation in which the state of the system is assumed to be coherent. With this 
assumption, each normal-ordered correlation function can be readily substituted by the 
corresponding product of field expectation values. The EOM read 



where s = (c s ) and a = (c a ) for brevity. The semiclassical approximation, in general, is 
valid for large photon numbers but it offers a satisfactory explanation of the dynamics 
of the system for smaller occupation as well, as we verify below with the numerical 
integration of the EOM (1261). 

The exact equations f |26l ) become nonlinear c-number equations for the expectation 
values s and a. Nonlinear equations may have multiple, stable and unstable steady-state 
solutions. A paradigmatic example of this behavior are the classical equations of the 
laser [62], where a solution which lacks coherent emission is possible at any pumping 
strength, but is unstable above a certain threshold. Here we find a comparable result, 



dt(cs) = [-i(6 d - 2J) - T](c s ) + n(c s c\c a ), 
d t (c a ) = [-i&d - r](c a ) -iQ- k{(c\c 8 + l)c a ) 



(26) 



d t s — i(Sd — 2J)s — (r — ft|a| 2 )s, 
d t a = —iS^a — id — [T + ft(|s| 2 + l)]a, 



(27) 
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Figure 4. Results (|3T|) of the semiclassical approximation to the EOM ([26]) . a) 
Population fraction n s /n in the symmetric mode as a function of cavity pumping Q and 
dissipation strength n. b) Population in the symmetric (solid line) and antisymmetric 
(dashed line) mode for n = V as a function of the pumping strength Q. The threshold 
(1281) is Q crit = 2.0 T. 



with a threshold value 

^ = ryi±S, (28) 

for the pumping strength Q. Below the threshold, the only stable steady-state solution 
to Eq. ( |27j) is (c.f. also |Appendix A[ ) 

|s|2 = »• |a|2 = <3 + (r + B y <29) 

and only the pumped antisymmetric mode is populated. For driving strengths above 
the threshold we obtain the stable solution 



,2_ ^ S l F , |2_ ^ 

"Vf^"^"^" 1, |a| ~ ^ + [r + /^(| 5 | 2 + i)] 2 ' (30) 

where the symmetric mode features nonvanishing occupation. Note that these equations 
depend on J only via the detuning 8^ which has the only effect to reduce the driving 
strength. Therefore, from now on we will restrict ourselves to 8^ = and J — >> 0. The 
stationary solution to the semiclassical EOM then simplifies to 

\s\ 2 = 0, \a\ 2 = ( r ^2 i if ^ < ^crit? ^gj^ 



crit • 



The behavior of the populations n 5 = |s| 2 and n a = |a| 2 is plotted in Fig. [Has a function 
of the pumping intensity and the effective decay k. Below threshold, the symmetric 
mode is empty and the antisymmetric mode occupation grows quadratically with the 
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pumping intensity. For Q > f2 cr i t , the antisymmetric mode population saturates, and 
the symmetric mode occupation grows linearly. Thus, if Q is sufficiently large, most of 
the photon population is transferred to the symmetric mode. This transition to a state 
in which photons 'macroscopically' occupy the symmetric superposition of two cavity 
modes is a direct consequence of the engineered nonlinear Liouvillian (|3]). 

The discussion of Eqs. ( 13T1) illustrates both, similarities and differences between 
the present setup and a laser (63J EH [65J 66j. In a single- or multi-mode laser, the 
light field is coupled to a pumped reservoir and the lasing threshold corresponds to the 
point where the total light intensity starts to grow. In our case, since the antisymmetric 
mode is driven directly by a classical field, the total photon population is always finite 
and grows monotonically. The threshold then marks the point where the many-body 
scattering mechanism contained in C K dominates, and a macroscopic transfer of photons 
from the antisymmetric to the symmetric mode occurs. In this sense, one can rather 
think of the antisymmetric mode as a photon reservoir from which a condensation of 
photons ('lasing') into the symmetric mode occurs. However, as we will see below, this 
analogy between lasing and condensation [TQl [IU [121 E3 EE] holds in our system only 
for certain parameter regimes, and in particular striking differences can be observed in 
the regime of low photon numbers. 

Finally, note that both for k as well as for k — >• oo the condensation is 
suppressed, and that the phase diagram in Fig. [4] is symmetric with respect to the ratio 
n/Y. For small k ) not enough photons are scattered into the symmetric mode, while for 
large k the antisymmetric mode is highly overdamped and the population of the whole 
system is small as in the first place. 

4-2. Exact diagonalization 

To support the results of the semiclassical approximation discussed above, we resort to 
the exact diagonalization of the system composed by two cavities. The time-evolution 
of the system for typical parameters features an initial transient in which the photonic 
population accumulates in the antisymmetric mode, which is directly pumped, and a 
later stage in which photons scatter into the symmetric mode. Finally, a steady state 
is reached in which a dynamical equilibrium takes place between the populations in the 
two modes. We remark that the population in the modes is continuously subjected to 
losses. 

The steady-state values of the populations n s = (cjc 5 ) and n a = (c\c a ) in the 
symmetric and antisymmetric mode, respectively, are shown in Fig. [5] as a function of 
the pumping strength Q and for two different values of k/Y. In the limit where 
we expect the transition to occur at large mode occupation numbers n a) we find that 
the exact solution matches qualitatively and quantitatively very well the semiclassical 
predictions. In particular, we see that above ^ crit the population n s grows linearly 
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Figure 5. Numerical solution to the exact equations of motion for L — 2 cavities. The 
occupation numbers n s and n a of the symmetric and antisymmetric mode are plotted as 
a function of the coherent driving strength Q for a) n/Y — 0.2 and b) k/Y = 5. In both 
plots #d = 0, and the dashed lines indicate the corresponding semiclassical predictions 
given by Eq. (|3T|) . The dotted lines show the two-photon correlation function of the 
symmetric mode = (n s (n s — l))/(n s ) 2 . 

with the driving strength, while the antisymmetric mode saturates at a value n a — Y/k. 
Differently from the prediction of the semiclassical analysis, however, n s does not vanish 
exactly below the threshold but it is smaller than n a . A study of the equal-time two- 
photon correlation g^ = (h s (h s — l))/ (n s ) 2 reveals that, when crossing fi cr i t , the photon 
statistics of the symmetric mode changes from a thermal (g^ w 2) to a Poissonian 
(g( 2 ) & 1) distribution. This aspect is in agreement with the standard lasing transition. 

In the opposite regime, k 3> T, Eq. f |3T1) predicts that the population n a of the 
driven mode is always less than one, and we do not expect the semiclassical analysis to 
give accurate results. Indeed, Fig. [5b shows that in this limit the transition is completely 
washed out and n s exceeds n a for all pumping strengths. Nevertheless, above ^ cr i t we 
still observe the characteristic saturation of n a and - apart from an additional offset - 
the correct linear scaling of the population in the symmetric mode. However, in strong 
contrast to the semiclassical regime, we find that the photon statistics of the symmetric 
mode is sub-Poissonian (g^ < 1) at low driving strengths and approaches the classical 
limit g^ = 1 from below. This anti-bunching effect can be understood from the fact 
that the effective damping of the c a mode is given by Y + K,((h s ) + 1). Therefore, the 
scattering of the first photon into the symmetric mode changes the damping significantly 
and suppresses the following repopulation of the 'reservoir mode' c a . The limit k 3> T 
and Q « Q cr it would then result in a 'single photon condensate' with n a <C n s ~ 1 and 
a non-classical statistics with g^ — >> 0.5. 
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4.3. Effective temperature 

Finally, let us briefly comment on the effect of the reverse photon scattering process 
~ k' in Eq. (ED- As we have discussed in Sec. [2J in the case of two isolated cavities the 
presence of the reverse photon scattering term can be interpreted as a finite temperature 
effect, and we can set k! = ^e _2J / feBTeff . By including this term, the EOM are 



These equations are of the same form of f |26l ). with the replacements r T + k! and 
k k — k! . Therefore, in the semiclassical regime we obtain the same physics, but with 
a renormalized pumping threshold 



In particular, for any fixed driving strength we can use Eq. (|33l) to define a critical 
effective temperature above which photon condensation does not occur. 

5. Photon condensation in a cavity array 

We turn now to the solution of Eq. ([2) in a lattice with L> 1 cavities. The central goal 
of this section is to show that the many-body Liouvillian (|3]) induces condensation 
of photons into the homogeneous lattice mode with zero momentum q. Following 
the scheme introduced in the previous section in the case of L = 2 cavities only, we 
consider the scenario in which each cavity mode q in the array, with £ = 1, . . .L, is 
coherently driven by a classical field with staggered amplitude Q# = (— In this way, 
the coherent driving acts on the edge of the Brillouin zone, while dissipation-induced 
condensation is expected to take place at the center. 

For simplicity, we assume periodic boundary conditions for the lattice. We rewrite 
Eq. (Ej) in terms of the annihilation operators c q = e iqi C£ of the photonic modes 

in momentum space, where q takes values — tt(L — 2)/L, . . . , tt(L — 2)/L, tt (given here 
in units of the inverse lattice spacing) in the discretized Brillouin zone. In the rotating 
frame with respect to the coherent driving, the Hamiltonian is then given by 



dt(c 8 ) 
dt(c a ) 



(r + K')(6a) -itl-(K- K')((C% + l)c a ). 



(32) 




(33) 



n c + nn = 5>d - 2 J(l + cos(q))]clc q + ^(^cj + fi^), (34) 



q q 

where Q q = e iq£ fli and 5^ = uj c + 2 J — uj^ is the detuning of the q = tt mode from 

the driving frequency. The LiouviUe operators are 




Photon condensation in circuit QED by engineered dissipation 



16 



and 

CkP= ^2 + * - <?2 - ^4)^(^1,^2,^3,^4) 



<?1<?2<?3<J4 
St 

In the latter equation, the function K plays the role of a scattering kernel, which has 



x [2clc ?2 pcJ Cq 4 c g4 cj c ?2 p c,j 4 cj c ?2 ]. (36) 



the general form 



4 

K(qi,q2,q3,q4) = j 



<?1 . <?2 . ?3 <?4 

k; cos — sm — sm — cos — 
2 2 2 2 

. / • Qi Q2 qs . <?4 
+ k sm — cos — cos — sm — 
2 2 2 2 J 



(37) 



In the following analysis, we will focus on the case k! = 0. Furthermore, as discussed 
above, the photon condensation effect does not depend on the tunneling amplitude and 
we can assume J — >> in the Hamiltonian. 

Solving the ME ([2j) in an extended array is in general a very demanding task, 
and exact diagonalization can only be provided for very limited system size (see our 
results in Fig. [6b for L = 4 cavities). For this reason, we resort in the following to 
two complementary approximations, which allow us to treat the full dynamics of the 
photonic population in the whole Brillouin zone. First we consider the equation of 
motion for the field expectation values within the semiclassical approximation, which 
was discussed in Sec. [4] for the case L = 2, and showed reliable results for large photon 
numbers. Then we treat the photonic population directly and derive a Boltzmann-like 
equation starting from the exact EOM for the two-point correlation functions. Within 
the latter approximation we are able to estimate the finite correlation length of the 
photonic field in the array. In this respect, this technique fills the gap between the 
semiclassical analysis of an infinitely extended array and the numerical solution to a 
small-size system. We emphasize that all techniques agree on the central results that 
we present here, i.e. the existence of a macroscopic photonic population at the center of 
the Brillouin zone. 



5.1. Semiclassical approximation 

For many cavities, the EOM for the cavity fields ip q = (c q ) (analogous to Eq. (1271)) read 
dt^q [-iSd ~ T - ac(1 - cosg)]^ itt q 

+ J 2^^i^«2^gi+ffl-g[ cos 9 ~ COS (?1 + ?2 - 9)], (38) 
Q1Q2 

where Q q = Sq^Q^. The numerical solution of this system of nonlinear coupled equations 
is straightforward. Fig. Efet shows the population n q = \ip q \ 2 on the discretized Brillouin 
zone for L = 22 cavities. We see that, in the initial stage of the dynamics, the coherently 
pumped mode at q = tt is populated on the time scale of T -1 . Later, photonic population 
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Figure 6. a) Numerical solution of the semiclassical equations of motion (|38j) . L = 22 
momentum states are used in the integration and the population n q = \ip q \ 2 is shown 
in the positive half q > of the Brillouin zone, b) One-particle density matrix (c q c q >) 
in momentum space, obtained with the numerical integration of the exact equations 
of motion for L = 4 cavities, integrated until t max T = 0.5. In both panels we use 
k = 5.0 T and fl q=7r = 3.5VLT. 



is transferred to the mode q = and the steady state is reached. A surprising feature is 
that the population of the intermediate modes in the Brillouin zone is almost vanishing 
(and indeed not visible in Fig. EH) throughout the time-evolution. This contrasts sharply 
with typical relaxation phenomena in the Brillouin zone (see e.g. Ref. [69J) where the 
excited population is expected to drift to lower energy states by interaction-induced 
relaxation. The reason for this behavior is the specific form of the engineered system- 
bath coupling C K . 

We can take advantage of the negligible population in the intermediate modes of 
the Brillouin zone and introduce a two-mode approximation, in which only the averages 
ipo and t/v of the modes q = and q = 7r, respectively, are assumed finite. The EOM 
(1381) then reduce to 

^o = -(r-4«:^|^ 7r | 2 )^o, m) 
^ 7r = -^ 7r -[r + 4K(i|^ | 2 + i/2)]^. 1 ; 

The structure of these equations is exactly the same as Eq. (1271) . with the |a| 2 /2 and 
\s\ 2 /2 replaced by the photon densities l^p/L and |^ | 2 /I/, respectively. The analytical 
discussion of these coupled equations follows the same path as that from Eq. f l27l ) to 
Eq. (ED). The threshold, above which photons start to condense into the q = mode, is 
^7r,crit = \/L/2 x \/r(r + 2k) 2 /2k. Note that compared to the results presented in the 
previous section, a factor 2ft, instead of ft, appears. This is simply due to the fact that 
in a ID array (with periodic boundary conditions) each cavity is dissipatively coupled 
to two neighboring sites. 
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In Fig. [6b we also present the results of the exact diagonalization of a small cluster 
of L = 4 cavities, for typical parameters comparable to the semiclassical results. We 
show the absolute values of the entries of the one-particle density matrix in momentum 
space, (c^ q c q f). The values along the diagonal are the populations n q) and the peak at 
q = clearly demonstrates substantial occupation of that mode. The non-diagonal 
entries, which represent the correlation between modes with different momentum, are 
small. These values cannot be obtained correctly from the semiclassical approximation, 
where correlations factorize and one would have |(cjc^/)| ~ y^V^V* ^ e a bsence 
of correlation between modes with different momentum can be understood, since the 
photonic scattering term originates from the non-unitary contribution to the ME (Ej). 
This fact suggests that an approximation scheme should hold, where the correlations 
between momentum modes are neglected from the beginning, and only the populations 
n q are the dynamical variables. The Boltzmann-like EOM that implements this scheme 
is presented below. 



5.2. Steady-state photon distribution 

The exact EOM for the populations n q = (c^ q c q ) are 

d t (c\c q ) = - 2[r + «(1 - cosq)]{c\c q ) + iQ*(c ? ) - ifi 9 (cj> 

+ T^ cos2 | sin2 | (6 ^- ) + lS((?) ' (40) 

where the term S(q) contains fourth-order correlation functions, which are not diagonal 
in momentum space, and reads 

<%) = j [ cos (^ ~ cos fai + q)] (c\ 2+qi _ q c\c qi c q2 ) + H.c. (41) 

<71<72 

One route to reduce the EOM to a manageable form is to implement a truncation in the 
infinite hierarchy of coupled correlation functions, expressing higher-order correlations 
in terms of products of lower-order ones. Here, we reduce the fourth-order correlation 
function to a product of second-order diagonal correlation functions according to the 
prescription 

(^l^L^S^l) — (^1,94^2,93 ~^ ^1,93^2,34) (^gi^<?l) (^2^2)' (42) 

This prescription is in fact a Hartree-Fock approximation [70], because we assume 
that the dominant contribution to the correlations arises in the density channel of the 
momentum modes. The motivations for this choice have been discussed above. Using 
this factorization we obtain the following simplified form for the scattering term 

-~l21 [ cos ^^ ~ cos (^)] (4i c 9i)( c K)- ( 43 ) 

<?1 
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Figure 7. Results of the numerical integration of the Boltzmann-like equations for the 
populations n q on a discretized Brillouin zone with L = 201 points, until t max T — 10.0, 
for k = 5.0 T, Q n = 3.5V^r (solid line in (a) and circles in (b)) and k = 1.0, 
Q n = 2.0\ /r LT (dotted line in (a) and crosses in (b)). Two Bose-Einstein equilibrium 
distributions are plotted in (b) (solid lines) for comparison. 



A self-consistent solution is to assume that (c q ) = for q ^ n. The EOM f |40l) become 
then a set of coupled, nonlinear, Boltzmann-like equations [69] for the populations 
only. However, on the edge of the Brillouin zone, the coherent pumping couples 
to the average field, which does not vanish. The numerical integration of the resulting 
equations is straightforward and produces the distribution n q in the whole Brillouin zone. 
A finer discretization in the Brillouin zone is possible, with respect to the semiclassical 
equations f l38l) . with a comparable numerical effort. Typical profiles of the steady-state 
distributions are shown in Fig. [71 Due to the increased resolution in momentum space, 
the broadening of the condensate peak at q = can be appreciated, and its dependence 
on the parameters is discussed below. 

To discuss the width of the central peak, which is related to the correlation length 
t; of the condensate (which we give here in units of the lattice spacing), we can restrict 
our analysis to the dynamics of the n mode and to the populations of modes with q w 0. 
For the coherence and population of the q = n mode we obtain 

d t (c n ) ~ -(i5 d + T + 2/c + 44)(0 " i^Tr, (44) 

and 

d t {clc n ) ~ -2 [r + 2k + 4«A/" ] (clc„) + i^(c n ) - iQ n (ct). (45) 
Here we have defined 

A/- = ij> g , (46) 
q 
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with Y2' q excluding q = 7r, and used the Hartree-Fock decoupling approximation, which 
for three operators reads ($ qi + q2 - q c qi c q2 ) ~ n qi 5(q 2 -q)(c q )+n q2 5(q 1 -q)(c q ). For 5 d = 0, 
the normalized steady-state population A/*^ = (cJ.c 7F )/L of the driven mode is then given 
by 

[r + 2«(1 + 2Af )f 
For the low-momentum modes we find 



K = ^^T7^ . (47) 



d t n g ~ - [2r - 8kA/; + k(1 + 2A/" + 2A4)g 2 ] n g + 8k(1 - <? 2 /4)AT 7r , (48) 
and the steady-state solution can be approximately written in the form 

< 49 > 

where 

8«JV; c2 _ k(1 + 2A/" + 2AQ 1 , . 

n °-2r-8^' e - 2r-s«A4 4- (50) 

Equation ( H9l) is a Lorentzian profile with broadening This scale plays the role 
of a correlation length since, in the translationally-invariant case, 

= £ ^e-^'{c\c q ,) c e-^'\c\c q ) oc e"l^. (51) 

Here we have used that, in momentum space, the single-particle density matrix is 
mainly occupied on the diagonal. This approximation is in line with the Hartree- 
Fock approximation used above and can be physically justified by the fact that low- 
momentum modes are populated by photons which are scattered incoherently from the 
q = TT mode. 

5.3. Discussion 

For the Lorentzian low-momentum distribution f l49l) we obtain Ao — ^o/£> an d we can 
interpret the quantity L x Ao as the total number of condensed photons, while the 
population L x A/*^ serves as finite photon reservoir. In steady state, the values of A/^, 
no and £ can be determined from the nonlinear Eqs. f |47l) and f |50l) . The stability of the 
low-momentum modes requires that A/*^ < Y /4ac, but for a driving strength above fi^crit 
the value of jV n will approach this limit. Therefore, in this regime 

where = |^7r|/^7r,crit > 1. We see that above the critical driving strength the number 
of photons in the low-momentum modes increases linearly with Cl^. At the same time 
the correlation length increases ^ f^, indicating an even stronger tendency for photons 
to occupy the q = mode. Note that while in the classical high-photon limit, f > ^, we 
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obtain the ratio Afo/Af n ~ (£\ — 1), in the opposite case, ^ > T, we get an additional 
enhancement factor, Ao/A/^ ~ (2^/r)(Q 7r — 1). Therefore, in this regime, very pure, 
dilute photon condensates with Ao ~ 1 and £ ~ (/^/4r)(^ 7F — 1)^ 3> 1 can be prepared. 
For typical parameters, the correlation length ranges from a few to a few tens of cavities, 
and can easily be tuned by adjusting the driving strength. This means that with realistic 
settings of around ten cavities, the variation of the correlation length could be studied 
in experiments. 

Finally, we find it worthwhile to compare the steady-state distribution of the non- 
equilibrium open system studied here to a thermal equilibrium gas of massive bosons 
with finite chemical potential. To do so, we compare the distribution f |49l ) with a 
Bose-Einstein distribution n^ E = l/(exp{(s q — ^)/fc B r} — 1) (c.f. Fig. [7b), where 
s q = — 2Jcos(g) in our case. Considering small momenta, a low q expansion of both 
distributions gives the following effective temperature and chemical potential for our 
photon condensate 

kBT 1 + n ° M = -2-I^lnfl + lY (53) 



j e ' j e v no, 

Note that the temperature T defined in this way is different from T e ^ associated to a 
finite ft', and arises purely from the competition between photon losses and equilibration 
processes in a driven system. Above threshold, where n ,^ > 1, we obtain the scaling 
k B T/J ~ 4/fi and k B T / J ~ 2T/(kCi) for the limits T > ft and T <C ft, respectively. 
This confirms the conclusion from above that, for ft/T —> oc, arbitrarily pure photonic 
condensates can be prepared with our scheme. 

As a final remark, we point out that the detuning 5^ of the driving field, which enters 
as an energy offset in the cavity Hamiltonian ( 1341) , affects the effective chemical potential 
ji only indirectly via a modification of f^crit- In particular, and in strong contrast to 
equilibrium systems, the sign of 5 a does not play a role in our non-equilibrium scenario 
and the chemical potential is related to the strength, rather than to the detuning of 
the driving field, as seen from the combined Eqs. f |53l ) and f |52l ). While for the present 
system this relation can be derived explicitly, it also suggests a general way to think 
about the effect of driving fields on the stationary states of open photonic quantum 
many-body systems. 



6. Conclusions and outlook 



In conclusion, we have analyzed a scheme to achieve condensation of microwave photons 
in an open and driven array of superconducting cavities. In particular, we have 
shown how the coupling to superconducting qubits can be used to engineer non-local 
and number-conserving dissipation processes for microwave photons, which mimic the 
coupling to an effective low-temperature bath. Under state-of-the-art experimental 
conditions, these processes can exceed the intrinsic losses in the system and produce 
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stationary many-body states of photons, which are independent of the details of the 
external driving fields. We have proposed a basic experimental setup for demonstrating 
a stationary condensate of photons in two or multiple coupled cavities and compared the 
properties of such an out-of-equilibrium photon condensate to an equilibrium scenario 
for massive bosons. 

More generally, our work illustrates the potential of superconducting cavity arrays 
for simulating various dissipative and out-of-equilibrium quantum many-body problems. 
The present example of photon condensation already shows how in engineered and 
fully controlled systems, an appropriately designed coupling to a bath can lead to non- 
trivial stationary states of an otherwise non-interacting photonic system, and could 
provide new insights for related experiments in condensed-matter systems. Beyond this 
basic example, the proposed scheme for engineering dissipation can be easily extended 
and combined with strong coherent photon-photon interactions. Compared to original 
ideas described in the context of cold atoms, superconducting circuits represent a 
complementary platform with, in many respects, additional flexibility to design non- 
standard dissipation processes, as well as to study open, non-equilibrium quantum 
many-body systems. 

Acknowledgments 

The authors would like to thank P. Zoller for proposing the initial direction of this 
project and further insightful discussions. We also thank F. Marquardt for valuable 
feedback on this work. This work was supported by the EU network AQUTE and the 
Austrian Science Fund (FWF) through SFB FOQUS and the START grants Y 591-N16 
(P.R.) and Y581-N16 (S.D.). 

Appendix A. Supplement to Section [4] 

In this Appendix we extend Section 0J Specifically, we study the case of incoherent 
pumping of the cavities and analyze the stability of the steady-state solutions found in 
that section. 

Incoherent pumping. It is important to notice that similar condensation effects to 
those shown in Section [4] can be achieved under incoherent pumping of the cavity modes. 
Neglecting the term Q (c\ + c a ) in the model, and including 



-pump 



(A.l) 



we arrive at the EOM 



d t (n s ) = 2p((h s ) + 1) - 2r(n s ) + 2/t(n a (n s + 1)). 
dt(h & ) = 2p«n a > + 1) - 2r(n a ) - 2/t(n a (n s + 1)) 



(A.2) 
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Figure Al. a) Fraction of photons in the symmetric mode, n s /n = (h s )/n as a 
function of the dissipative rate ft, for incoherent pumping p/Y = 0.5. b) Threshold 
given by Eq. (|28j) . below which ([29]) is the only stable solution and above which ([30]) 
is the only stable solution. 



Solving for the steady-state solution of these equations within the semiclassical 
approximation (h s h a ) & (h s )(h a ), we find the condition ^ = being n = (h s ) + (n a ) 
the total number of particles in the steady state. Substituting then (h a ) = n — (n s ) = 

— (n s ) in the first steady-state equation, for p ^ T, we find a smooth crossover 
between an equal population of symmetric and antisymmetric modes for n <C T, and 
condensation to the symmetric mode for k 3> p, T (see Fig. lAlb ). Interestingly, for 
p —> Y, however, all photons occupy the symmetric mode for all values of k ^ 0. In the 
latter case, heating (incoherent pumping) compensates cooling, and only the engineered 
dissipative mechanism survives. 

Stability of the steady-state solutions. The stability of the different solutions can 
be checked by linearization of the EOM around the steady-state values f |29l ) and f |30l ). 
To this end we linearize the EOM (1271) as s — >■ <s + 5s, a — » a + Sa (similarly for 
the conjugate fields), being s , &o the steady-state solutions and 5s, 5a representing 
fluctuations. This gives ^(5s, 5a, 5s*5a*)T = J (5s, 5a, 5s*5a*)T with the 4x4 Jacobian 
matrix J: 

( — (r — ft|a | 2 ) ^l^ol^o ft|so|a \ 

—tv\so\ao i5 c — [r + ft(|<s | 2 + 1)] ~~ ^l^ol^o 
ft|so| a o — (r _ ^| a o| 2 ) ^|^o| a o 

y -K\s \aQ -K\s \aQ -i5 c - [Y + ft(|s | 2 )] J 

whose eigenvalues will have negative real part in the region of stability. Evaluating 
the eigenvalues of this matrix at the steady state solutions, we find that the steady- 
state solutions presented in Section [4] are stable in the whole range of parameters. In 



Photon condensation in circuit QED by engineered dissipation 



24 



Fig. lAlb we show the threshold f2 crit given by Eq. (1251) . below which ( |29l) is the only 
stable solution and above which f |30l) is the only stable solution. 
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